# 行为数据线性混合效应模型 - Homograph子模型分析（使用词频比例）
# Behavioral Data Linear Mixed Effects Model - Homograph Submodel Analysis (using frequency proportions)

# Clean workspace
rm(list = ls())  # Remove all objects
graphics.off()    # Close all graphics devices
cat("\014")      # Clear console

# Install required packages if not already installed
if (!require(tidyverse)) install.packages("tidyverse")
if (!require(readxl)) install.packages("readxl")
if (!require(writexl)) install.packages("writexl")
if (!require(moments)) install.packages("moments")
if (!require(nortest)) install.packages("nortest")
if (!require(lme4)) install.packages("lme4")
if (!require(car)) install.packages("car")
if (!require(lattice)) install.packages("lattice")
if (!require(ggplot2)) install.packages("ggplot2")
if (!require(DHARMa)) install.packages("DHARMa")
if (!require(emmeans)) install.packages("emmeans")
if (!require(sjPlot)) install.packages("sjPlot")

# Load required libraries
library(tidyverse)  # for data manipulation
library(readxl)     # for reading Excel files
library(writexl)    # for writing Excel files
library(moments)    # for skewness and kurtosis calculations
library(nortest)    # for Anderson-Darling normality test
library(lme4)       # for mixed effects models
library(car)        # for additional diagnostics
library(lattice)    # for plotting
library(ggplot2)    # for visualization
library(DHARMa)     # for residual diagnostics
library(emmeans)    # for simple effects analysis
library(sjPlot)     # for model visualization

# Function to load data
load_behavior_data <- function(file_path, sheet_number = 3) {  # Default to load the 3rd sheet
  # Determine file extension
  file_ext <- tools::file_ext(file_path)
  
  # Load data based on file type
  data <- if (file_ext == "csv") {
    read.csv(file_path)
  } else if (file_ext %in% c("xlsx", "xls")) {
    read_excel(file_path, sheet = sheet_number)  # Specify sheet
  } else {
    stop("Unsupported file format. Please provide a CSV or Excel file.")
  }
  
  # Convert to tibble for better handling
  data <- as_tibble(data)
  
  # 处理变量名：将空格替换为下划线
  names(data) <- gsub(" ", "_", names(data))
  
  # Ensure proper column types
  data <- data %>%
    mutate(
      Subject = as.factor(Subject),
      Semantic_Relatedness = as.factor(Semantic_Relatedness),
      Priming_Word_Type = as.factor(Priming_Word_Type),
      Group = as.factor(Group)
    )
  
  return(data)
}

# Set input file paths
behavior_file <- "source_data/Ch_Ja_Lex_Behavior_2Group_29Sub_for R.xlsx"
freq_file <- "results/word_frequency/word_frequency_processed.xlsx"

# Create result folder name from input file name
result_folder <- tools::file_path_sans_ext(basename(behavior_file))
result_folder <- gsub("_for R", "", result_folder)  # Remove "_for R" suffix
result_path <- file.path("results", result_folder, "mixed_effects_homograph_prop")

# Create result folder
if (!dir.exists(result_path)) {
  dir.create(result_path, recursive = TRUE)
  print(paste("Created result folder:", result_path))
}

# Load behavior data
print("Loading behavior data...")
behavior_data <- load_behavior_data(behavior_file, sheet_number = 3)
print("Behavior data dimensions:")
print(dim(behavior_data))
print("Number of homograph trials:")
print(sum(behavior_data$Priming_Word_Type == "homograph"))

# Load word frequency data
print("Loading word frequency data...")
word_freq_data <- read_excel(freq_file)
print("Word frequency data dimensions:")
print(dim(word_freq_data))
print("Word frequency data column names:")
print(names(word_freq_data))

# Filter homograph data and merge with word frequency
print("Processing data...")
homograph_data <- behavior_data %>%
  filter(Priming_Word_Type == "Homograph") %>%
  left_join(word_freq_data, by = c("Prime_word" = "Word"))  # 合并词频数据

print("Homograph data dimensions after merge:")
print(dim(homograph_data))
print("Number of NA values in key columns:")
print(colSums(is.na(homograph_data)))

# Make sure log_RT exists
if (!"log_RT" %in% names(homograph_data)) {
  print("Creating log-transformed RT...")
  homograph_data <- homograph_data %>%
    mutate(log_RT = log(RT))
}

# 创建iterm变量（组合Prime word和Target word）
homograph_data$iterm <- interaction(homograph_data$Prime_word, homograph_data$Target_word)
homograph_data$iterm <- as.factor(homograph_data$iterm)

# 检查模型变量是否存在
print("Checking model variables:")
print(names(homograph_data))
print("Checking for NA values in model variables:")
print(colSums(is.na(homograph_data[c("log_RT", "Semantic_Relatedness", "Group", "Subject", "iterm")])))

# 检查因子水平
print("Checking factor levels:")
print("Semantic_Relatedness levels:")
print(levels(homograph_data$Semantic_Relatedness))
print("Group levels:")
print(levels(homograph_data$Group))

# 设置模型公式
model_formula <- as.formula("log_RT ~ Semantic_Relatedness * Group + 
                          (1 + Semantic_Relatedness | Subject) + 
                          (1 + Group | iterm)")
model_formula_text <- "log_RT ~ Semantic_Relatedness * Group + 
                          (1 + Semantic_Relatedness | Subject) + 
                          (1 + Group | iterm)"

print(paste("Model formula:", model_formula_text))

print("Fitting linear mixed-effects model...")
mixed_model <- lmer(model_formula, data = homograph_data)
print("Model fitting completed")

# 模型诊断
print("Performing model diagnostics...")

# 1. 检查模型收敛
print("Model convergence check:")
print(isSingular(mixed_model))

# 2. 检查随机效应
print("Random effects:")
print(ranef(mixed_model))

# 3. 检查固定效应
print("Fixed effects:")
print(summary(mixed_model))

# 4. 残差诊断
print("Residual diagnostics:")
residuals <- resid(mixed_model)
print(summary(residuals))

# 5. 正态性检验
print("Normality test for residuals:")
print(ad.test(residuals))

# 6. 方差膨胀因子检查
print("Variance Inflation Factors:")
print(vif(mixed_model))

# 保存模型结果
print("Saving model results...")
model_summary <- summary(mixed_model)

# 保存模型摘要到txt文件
sink(file.path(result_path, "model_summary.txt"))
cat("线性混合效应模型结果分析报告\n")
cat("================================\n\n")

# 1. 模型概述
cat("1. 模型概述\n")
cat("   模型公式：", model_formula_text, "\n")
cat("   数据观测数：", nrow(homograph_data), "\n")
cat("   受试者数量：", length(unique(homograph_data$Subject)), "\n")
cat("   项目数量：", length(unique(homograph_data$iterm)), "\n\n")

# 2. 固定效应结果
cat("2. 固定效应结果\n")
cat("   以下报告所有固定效应的估计值、标准误、t值和p值：\n")
cat("   显著性标记说明：\n")
cat("   p < 0.1: ·\n")
cat("   p < 0.05: *\n")
cat("   p < 0.01: **\n")
cat("   p < 0.001: ***\n\n")

# 从模型摘要中提取固定效应
fixed_effects <- fixef(mixed_model)
se <- sqrt(diag(vcov(mixed_model)))
t_values <- fixed_effects / se
p_values <- 2 * (1 - pnorm(abs(t_values)))

# 添加显著性标记
significance <- sapply(p_values, function(p) {
  if (p < 0.001) "***"
  else if (p < 0.01) "**"
  else if (p < 0.05) "*"
  else if (p < 0.1) "·"
  else ""
})

results_df <- data.frame(
  Term = names(fixed_effects),
  Estimate = fixed_effects,
  SE = se,
  t_value = t_values,
  p_value = p_values,
  Significance = significance
)

print(results_df)
cat("\n")

# 3. 随机效应结果
cat("3. 随机效应结果\n")
cat("   报告随机效应的方差组分：\n\n")
print(VarCorr(mixed_model))
cat("\n")

# 5. 结果解释
cat("5. 结果解释\n")
cat("   5.1 主效应分析\n")
cat("   - Group主效应：", ifelse(results_df$p_value[results_df$Term == "Groupexperimental"] < 0.05, "显著", "不显著"), 
    results_df$Significance[results_df$Term == "Groupexperimental"], "\n")
cat("   - Semantic_Relatedness主效应：", ifelse(results_df$p_value[results_df$Term == "Semantic_RelatednessSemantic Unrelated"] < 0.05, "显著", "不显著"),
    results_df$Significance[results_df$Term == "Semantic_RelatednessSemantic Unrelated"], "\n\n")

cat("   5.2 交互作用分析\n")
cat("   - Group × Semantic_Relatedness：", ifelse(results_df$p_value[results_df$Term == "Groupexperimental:Semantic_RelatednessSemantic Unrelated"] < 0.05, "显著", "不显著"),
    results_df$Significance[results_df$Term == "Groupexperimental:Semantic_RelatednessSemantic Unrelated"], "\n\n")

# 6. 随机效应解释
cat("6. 随机效应解释\n")
cat("   - Subject随机截距和Semantic_Relatedness随机斜率的方差表示个体间的变异程度\n")
cat("   - Item随机截距和Group随机斜率的方差表示项目间的变异程度和Group效应在项目间的变化\n\n")

# 7. 总结
cat("7. 总结\n")
cat("   根据以上分析结果，主要发现：\n")
cat("   - Group效应", ifelse(results_df$p_value[results_df$Term == "Groupexperimental"] < 0.05, "显著", "不显著"),
    results_df$Significance[results_df$Term == "Groupexperimental"], "\n")
cat("   - Semantic_Relatedness效应", ifelse(results_df$p_value[results_df$Term == "Semantic_RelatednessSemantic Unrelated"] < 0.05, "显著", "不显著"),
    results_df$Significance[results_df$Term == "Semantic_RelatednessSemantic Unrelated"], "\n")
cat("   - Group × Semantic_Relatedness交互作用", ifelse(results_df$p_value[results_df$Term == "Groupexperimental:Semantic_RelatednessSemantic Unrelated"] < 0.05, "显著", "不显著"),
    results_df$Significance[results_df$Term == "Groupexperimental:Semantic_RelatednessSemantic Unrelated"], "\n\n")

cat("注：所有效应的具体方向和大小请参考上述详细统计结果。\n")
cat("显著性标记说明：\n")
cat("   · 表示 p < 0.1（边缘显著）\n")
cat("   * 表示 p < 0.05（显著）\n")
cat("   ** 表示 p < 0.01（非常显著）\n")
cat("   *** 表示 p < 0.001（极其显著）\n\n")
cat("分析日期：", format(Sys.Date(), "%Y年%m月%d日"), "\n")
sink()

# 保存结果到Excel
write_xlsx(results_df, file.path(result_path, "model_results.xlsx"))

# 简单效应分析
print("Performing simple effects analysis...")

# 1. Group在Semantic_Relatedness各水平上的简单效应
group_simple_effects <- emmeans(mixed_model, ~ Group | Semantic_Relatedness)
group_simple_effects_df <- as.data.frame(group_simple_effects)
write_xlsx(group_simple_effects_df, file.path(result_path, "group_simple_effects.xlsx"))

# 2. Semantic_Relatedness在Group各水平上的简单效应
semantic_simple_effects <- emmeans(mixed_model, ~ Semantic_Relatedness | Group)
semantic_simple_effects_df <- as.data.frame(semantic_simple_effects)
write_xlsx(semantic_simple_effects_df, file.path(result_path, "semantic_simple_effects.xlsx"))

# 3. 成对比较
pairwise_comparisons <- pairs(emmeans(mixed_model, ~ Group * Semantic_Relatedness))
pairwise_comparisons_df <- as.data.frame(pairwise_comparisons)
write_xlsx(pairwise_comparisons_df, file.path(result_path, "pairwise_comparisons.xlsx"))

# 将简单效应分析结果添加到模型摘要中
sink(file.path(result_path, "model_summary.txt"), append = TRUE)
cat("\n8. 简单效应分析\n")
cat("   8.1 Group在Semantic_Relatedness各水平上的简单效应\n")
print(group_simple_effects)
cat("\n")

cat("   8.2 Semantic_Relatedness在Group各水平上的简单效应\n")
print(semantic_simple_effects)
cat("\n")

cat("   8.3 成对比较结果\n")
print(pairwise_comparisons)
cat("\n")

cat("注：简单效应分析用于进一步解释Group和Semantic_Relatedness的交互作用。\n")
cat("   具体结果请参考生成的Excel文件。\n\n")
sink()

# 可视化
print("Creating diagnostic plots...")

# 1. 残差图
png(file.path(result_path, "residual_plot.png"))
plot(mixed_model)
dev.off()

# 2. QQ图
png(file.path(result_path, "qq_plot.png"))
qqnorm(residuals)
qqline(residuals)
dev.off()

# 3. 随机效应图
png(file.path(result_path, "random_effects_plot.png"))
dotplot(ranef(mixed_model))
dev.off()

# 4. 固定效应图
png(file.path(result_path, "fixed_effects_plot.png"))
plot_model(mixed_model, type = "pred", terms = c("Semantic_Relatedness", "Group"))
dev.off()

print("Analysis completed successfully!") 